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Abstract 

This work addresses the inverse problem of electrocardiography from a new perspective, 
by combining electrical and mechanical measurements. Our strategy relies on the defini¬ 
tion of a model of the electromechanical contraction which is registered on EGG data but 
also on measured mechanical displacements of the heart tissue typically extracted from 
medical images. In this respect, we establish in this work the convergence of a sequential 
estimator which combines for such coupled problems various state of the art sequential 
data assimilation methods in a unified consistent and efficient framework. Indeed we ag¬ 
gregate a Luenberger observer for the mechanical state and a Reduced Order Unscented 
Kalman Filter applied on the parameters to be identified and a POD projection of the 
electrical state. Then using synthetic data we show the benefits of our approach for the 
estimation of the electrical state of the ventricles along the heart beat compared with 
more classical strategies which only consider an electrophysiological model with EGG 
measurements. Our numerical results actually show that the mechanical measurements 
improve the identifiability of the electrical problem allowing to reconstruct the electrical 
state of the coupled system more precisely. Therefore, this work is intended to be a first 
proof of concept, with theoretical justifications and numerical investigations, of the ad¬ 
vantage of using available multi-modal observations for the estimation and identification 
of an electromechanical model of the heart. 


1. Introduction 

In the last few years, more and more attention has been paid to the problem of state 
and parameters identification for complex three-dimensional models used in biomedical 
applications. Several works can be cited: for example in cardiac electrophysiology mi 
m [13 EZ] or in cardiac mechanics misiiisa, or also in hemodynamics 13 nnum 
E3- The observations available in these contexts are often multiphysics since several 
modalities can be used simultaneously: electrocardiograms, electrograms, MRI, GT scan, 
flow measurements with ultrasound, pressure measurement with catheters, myocardium 
thickness measurement with piezoelectric sensors, etc. Up to now, the multiphysics 
nature of these problems has sometimes been taken into account in the direct models 
but it has rarely been used in the inverse problems. 

Preprint submitted to Elsevier March 24, 2015 




The first purpose of the present study is to exhibit an example where multiphysics 
observations actually improve the identihability of a coupled system. More precisely, 
considering an electromechanical model of the heart, it is shown that the estimation 
of an electrical parameter is improved if the electrical observation, namely the elec¬ 
trocardiogram, is enriched with mechanical observations, namely the movement of the 
myocardium. In this model, the electrophysiology acts as an input for the mechanics, but 
the electromechanical feedback is neglected m- The coupling is therefore only one-way. 
The second purpose of this article is to show that a reduced order filtering strategy is 
well-suited to this class of multiphysics problems. Optimal filtering, like Kalman filter 
and its nonlinear extensions, is known to be very efficient but also too expensive to be 
used for large problems like those considered here. An effective strategy consists in using 
this kind of methods for the parameters only, and to address the uncertainties on the 
state variables through a less expensive approach. This strategy has been successfully 
used for example in mechanics where the state variable was handled with a Luenberger 
filter [33] ■ In this paper, a similar strategy is proposed to filter the electrical state vari¬ 
ables, by exploiting the special structure of the one-way coupled problem. But due to the 
structure of the equations, a Luenberger approach is not straightforward in electrophysi¬ 
ology. A different reduced order method, based on Proper Orthogonal Decomposition is 
then proposed, as it was done for cardiac mechanics in m- 


1.1. Background and related work 

The goal of the inverse problem of electrocardiography, also called cardiac electrical 
imaging, is to reconstruct the electrical activity of the heart from body surface potential 
maps. Various strategies have been proposed since four decades. All of them assume that 
measurements of the electrical potential ux are available on parts of the torso boundary 

The different strategies can be distinguished by the cardiac electrical source models 
they rely on. One of the first approaches was to estimate equivalent electrical dipoles 
[25l l89] . Another popular approach is to estimate the heart surface potential, usually 
called epicardial potential (even though pericardial potential would be more appropriate 
as noted in |3H|)- The potential mt within the torso is assumed to be solution of the 
Poisson problem: 

fV-(^^-VMT) = 0, inD- 
= Ue, on dflg 

where dfig denotes the boundary of the heart and Dj. is the electrical conductivity of 
the torso. The inverse problem then consists in estimating Ue on dflg (see e.g. |H 
da I61jb This problem being notoriously ill-posed, various regularizations have been 
proposed: Tikhonov [37], the use of temporal information |211llH|i truncated Singular 
Value Decomposition or truncated Total Least Square m- 

Another approach to address the inverse problem of electrocardiography is to consider 
the following equation within the heart Oq: 


—V • ((D. + D^) ■ VMe) = Y • (Dj • Vum), in Dq, 


which is one of two equations of the bidomain model (see Section 2.1). Then, instead 
of estimating the epicardial potential on the surface dflg, the goal is to estimate the 
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transmembrane potential Vm within the heart. In [38) . the two approaches were com¬ 
pared: investigating the null-space of the inverse problem, the authors concluded that the 
transmembrane potential-based formulation is more promising because it is based on a 
stronger biophysical a priori. As the epicardial potential approach, the transmembrane 
potential approach is ill-posed and must be regularized. In [34] . four variants of Lp'- 
Tikhonov regularization are compared. In [45] . i/^-Tikhonov regularization is used with 
a prior taking two different homogeneous values in the myocardium depending on cardiac 
phase (plateau or rest values). The inverse problem, formulated in a PDE-constrained 
framework, is addressed by directly solving the optimality saddle-point problem. In |46j . 
the estimation of the transmembrane potential is combined with a level set technique 
to efficiently identify the location of a myocardial infarction. In |62j . the approach of 
|45| is generalized to more general objective functions and constraints in order to iden¬ 
tify ischemic regions (characterized by lower amplitude during the plateau phase). Two 
different regularizations are investigated: the Tikhonov regularization, which is found to 
overestimate ischemic regions but with good sensitivity, and the Total Variation regular¬ 
ization which is found to underestimate ischemic regions but with high specificity. 

In [54] . the authors note that the usual regularization techniques have no physical 
ground. Instead, they propose to regularize the inverse solution with the monodomain 
equations coupled to the Fenton-Karma ionic model. The strategy proposed in the 
present paper has some similarities with this approach. We rely on a full electrophysio- 
logical model of the action potential coupled to the Poisson problem Q to estimate the 
solution of the heart electrical activation. This physical model is personalized on the fly 
with respect to its parameters in order to adapt it to a specific patient. Furthermore, 
we propose an additional step of modeling by considering the mechanical response to 
the electrophysiological activation, so we are able to also integrate mechanical measure¬ 
ments. Indeed, we believe that multimodal observations improve the identifiability of 
the complete model and therefore improve the quality of the electrical and mechanical 
state reconstruction. Another originality of our work is the use of a sequential data as¬ 
similation strategy that is adapted to a coupled electromechanical evolution model. Here 
we demonstrate how state-of-the-art gain filter on the electrophysiological model and on 
the mechanical model can be aggregated to propose a joint gain filter for the coupled 
problem. 

1.2. Organization of the present work 

The paper is organized as follows. In Section the electromechanical model is pre¬ 
sented. In Section]^ the observation operators - namely the measurements - are detailed 
for the electrical and the mechanical variables. In Section the general notions of data 
assimilation, optimal filtering and reduced order filtering are reviewed. Although the 
algorithms of this section are not new, their presentation differs from what is most often 
done in the literature since a purely deterministic description is adopted. In Section 
the algorithms used for the electromechanical problem are proposed and analyzed. In 
Section numerical experiments based on synthetic data are presented. The main pur¬ 
pose is to estimate a non-homogeneous parameter of the electrical model using electrical 
and mechanical observations. 
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2. Models 


We present the models in a time and space continuous context before entering into 
the discretization details and their numerical implementation. Concerning the continuous 
context, we denote by an underline character any vector field of and two underlines 
any second-order tensor. 


2.1. Electrophysiology 

A widely accepted model of the macroscopic electrical activity of the heart is the so- 
called bidomain model [HI |56l |59| [60] • It consists of two degenerate parabolic reaction- 
diffusion PDEs which describe the dynamics of the averaged intra- and extracellular 
potentials Mi and Me, coupled to a system of ODEs defining an ionic model. This model 
is related to the chemical dynamics of the myocardium cell membrane, in terms of the 
(vector or scalar) variable w representing the distributed ion concentrations and gating 
states, or a phenomenological counterpart. The model reads 

T '^)J ^ ' (^.Vm,) = Urn-Iapp; ’ 

C^micm^t'Cm T T — Clm^appi ill Hq ’ 

dtw + g{vm,w) = 0, in flo, 


where = Ui — represents the transmembrane potential, is the membrane capac¬ 
itance per unit area, Om is a constant representing the rate of membrane area per unit 
volume, D_.,D are the intra- and extra-cellular conductivity tensors, Japp is an external 
volume current. In ([^ the function g represents an ionic model. In this article, the 
Mitchell-Schaeffer ionic model [40] is considered, with the same rescaling as in [10]. It is 
a reduced complexity model capable of integrating relevant phenomenological properties 
of the ventricle cell membrane: 
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where Egate, Tin, Tout, Topen, Tdose are given constants and Enin and Enax are scaling 
constants (typically Emin = —SOrnE and Emax = 20mE). 

On the boundary dflg, we have E. ■ Vui ■ n = 0, and the heart is assumed to be 
isolated, • Vue ■ n = 0, as often done in the literature [HES]. For well-posedness, the 
condition J^h Me = 0 is enforced. 
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Summing and subtracting the first two equations of the system reads 


Y ' (—[ ' Y (^e T ^m)) — ^m4pp; 

in fig, 

—Y • ((Yj + Yg) ■ Y^e) — Y • (Yj ■ Y^^m) = 0, 

in fig. 

dtw + g{vm,w) = 0, 

in fig. 

(iZg • YUe) • Zt = 0, 

on 90 q, 

(4; • Y (We + Vm)) • n = 0, 

on 90 q. 


( 4 ) 


A Pi finite element discretization of the potential variables is used - see the corre¬ 
sponding refined electrical mesh in Figure]^ This leads to a discretization space with 
A^® degrees of freedom (i.e. V% ~ ). To each field - for instance Vm - is associated its 

corresponding approximation, denoted with an index h — for example Vmh- Equivalently, 
this approximation can be represented by its corresponding vector of degrees of freedom 
written with a vector in uppercase letter - i.e. Vm- The finite element vector of degrees 
of freedom or the discretization of linear form are then defined with vectors in straight 
uppercase letter whereas the finite element matrix operator with a bold uppercase letter. 
With this notation: 


V17^ e V^, 




V. 


mh^h 








dVl, 


VC/^ e VI, U'^VitVm = [ R - YVmh ■ YuJ, 

Jno 

e ^ (4pp - lionivmh, Wh)) < dn ~ (tpp - W)) 

and if a spatial discretization of the internal variable is done by node 

yU^GVl f giv,nh,Wh)uUnc^U^M^G{V^,W) 

where the applied current and the ionic variables have been interpolated, allowing to 
define the ionic variables at the nodes instead of the quadrature points. Finally, after 
spatial discretization system 0 reads 

+ K®(f4 + Ue) = Mf (4pp - Iio„(4n, W)), 

(Kf+ K|)C/e + K®t4 = 0, (5) 

W + G{V^,W) = 0. 

The state of this system is A® = ^ ^ ^ while Ue appears as an auxiliary variable verifying 
the static equilibrium ([^2 with Vm- This can be summarized by 

'a® = A®(A®,6»®), 

A®(0) = A® + Cxe, (6) 

= + 

where 0® denotes the vector of parameters of the electrophysiology model, and Cx® 
the uncertainty on the parameters and the initial condition respectively. 
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2.2. Mechanics 

The heart domain is denoted by n^{t) at any time t. This domain is the image of a 
reference conhguration fig through the solid deformation mapping cj) 

n^x[o,T]^n^{t), 

^ X = ^ + y(^, t) 

where y is the solid displacement. The solid velocity is given by w = y. The deformation 
gradient F is given by F{^,t) = = 1 + ^y, and its determinant is denoted by 

J. The right Cauchy-Green deformation tensor is defined by C = ■ F, the Green- 

Lagrange tensor by e = ^ (G — 1) = ^ (V^y + (+ (Y^y)^ ■ Y^y), and its linearization 

by 1 = W^^y + (.^^yy)- 

The mass per unit volume is denoted by p and the Cauchy stress tensor by tr. In the 
reference configuration, the first and second Piola stress tensor are respectively defined 
by T = Jq_ ■ F_~'^ and S = F_~^ ■ T_ = JF~^ ■ a ■ The constitutive law is assumed 

to be a combination of a hyperelastic law of potential W, a viscous component chosen 
proportional to the strain rate e, and an active part along the hber direction r represented 
by 3 internal variables which are Cc the active strain, kc the active stiffness and Tc the 
associated active stress m- 

dW 

Yfe Sc, fee, Tc) = -^ig) + Vsg + cr^oiec, kc, Tc)t O r, (7) 


with (Tin = (F: + MGc), where these 3 internal variables rely on a chemically- 

controlled constitutive law describing the myofibre mechanics [3111 [ 17 ]: 

f dtkc =-i\u\+a\ec\)kc + noko\u\^ in fig 

\ ( 8 ) 

[ dtTc = -(|m| -I- a \ec\)Tc + Cckc + noCTo |u|+ in fig 

with a, fco, (To given parameters, ng a function of Cc accounting for the Frank-Starling 
effect and u directly related to the electrical activity of the heart by 


u{t) = avm{t) -f h 


where a and b are two scaling parameters. 

Goncerning the boundary conditions, the external organs are modeled by visco-elastic 
boundary conditions on a sub-part of the epicardium: Tg-n = ksy + CsP on r„(t). A uni¬ 
form pressure is enforced on the left and right endocardium: = Pv,iilt on r„_i(t), i = 

{1,2}. In summary, the mechanical problem reads 


dty = 

-- V, 

in fig 

pdtv 

-Y-(Y) = 0, 

in fig 

dike ■■ 

= -( u -1- a\ec\)kc + nofco u ^ , 

in fig 

dtTc -- 

= -(|u| -1- a \ec\)Tc + Cckc + noag |m|_^ 

, in fig 

T-n 

= key + CsV, 

on r„ 

T-n 

— JVi, jF * n. 

on Teg 

T-n 

= 0, 

on i9fl| 


( 9 ) 


on i9fl5\((Uirc,i) ur„) 


6 







with the constitutive law (Hit- 
Remark 1 

Several improvements can be formulated on this model: other active or passive constitu¬ 
tive laws, more sophisticated boundary conditions for the tethering of the myocardium 
or hemodynamics, etc. The estimation strategy presented here can be adapted to any 
of these improvements since we only refer to state-space model description to formulate 
our data assimilation methods. 

The system is discretized with a Pi finite element - see the corresponding mechanical 
mesh in Figure [^- with a 5% compressibility acceptance in order to avoid any numerical 
locking. Concerning the fibre directions, we prescribed them on each point of the mesh 
with an elevation angle varying from -60 degrees to 60 degrees through the myocardium 
thickness. The discrete system is based on the variational formulation associated with 


€ V“, j pdtv ■ 2 ^ dV. + j E(e, e, Cc, A:c, Tc) : dye • u*'dn 

J sIq j 

-b f (ksU + Cgv) ■ dT =j Jpy{F~'^ ■ n) ■ dr (10) 
dr„ - , Jr.,. - 

Using the same convention as for the electrophysiological model discretization, the mass 
operator is defined by 


VU^U‘’^M“U= / pdtVf,-vidn, 


/a? 


the stress residual by 

^V\ V^^K^{Y,V)= f 

some weighted mass operators on the boundary by 


and a following pressure operator by 

■n)-vldr. 

r n . i 


The internal variables ec,kc,Tc are gathered in a vector which is discretized at the 
integration points. A vector of degrees of freedom ic is associated with the discontinuous 
field (ec,/i, kc^h, Tc^h) and the model Q is discretized by ic = B“(?c, Y, V). The complete 
mechanical model is spatially discretized into 


[f = U, 

< + K“(y, V, I,) + M“ (H) 

I t; = B“(z;,y,u). 
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The state of this system is X™ = T j • Denoting by 0™ the set of the parameters 

characterising the mechanics, affected by an a priori uncertainty and with the a 
priori uncertainty on the initial condition, the mechanical system reads 

X“(0)=X“ + Cx», (12) 

where the electrical variable X®, solution to can be seen as an input. 

2.3. Electromechanical coupling 

From a computational point of view, the electromechanical problem consists in solv¬ 
ing two coupled systems. We chose to keep the two sub-systems in independent solvers. 
This choice allows us to use legacy codes, to make their maintenance easier and to take 
advantage of the specific numerical methods adapted to each physical compartment. The 
coupling algorithm sketched in Figure is handled by a “master” code which exchanges 
the heart displacements and the transmembrane potential with the electrical and mecha¬ 
nical software. In this work, it is assumed that there is no electromechanical feedback. 
The transmembrane potential Vm is sent to the mechanical problem. Then, the one-way 
coupling is performed through the quantity u = avm + b which triggers the mechanical 
contraction via a change in the active stiffness kc and in the active stress Tc, see (|^. 



Figure 1: Master - Slave coupling 


Because of the characteristics of the action potential propagation, the space and time 
steps required by the electrophysiology are typically smaller than those needed by the 
mechanics. To ensure accuracy at a reasonable computational cost, each sub-problem is 
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solved with its own space and time step. The data are transmitted by the master code 
at some check-point in time (see Figure and interpolated in space from a mesh to 
another (see Figure^. Eventually, a typical complete direct simulation of our model is 
presented in Figure 



Figure 2: Mechanical mesh (Left, 4907 nodes, 18193 elements) and electrical mesh (center, 108112 nodes, 
541994 elements) and thorax mesh (right, 229782 nodes, 1250072 elements) 


The inverse problem will be addressed on the discrete formulations. The state and 
the parameters are defined as the combination of the electrical and the mechanical ones 
^ ~ ix^) ’ ^ ~ ie^) ■ Defining A(X, 6*) = ^ e")) ’ state variable follows 

the dynamics 


X = A{X,e) (13) 

The initial condition includes the quantities and which model the uncertainties 
on the initial condition and on the parameters respectively. When necessary, a trajectory 
will be denoted by (*<)]. This system can also be written in an augmented form by 
concatenated the state and parameters in = (^). This augmented state satisfies 
X = *’A(A). The real trajectory is assumed to be a solution of the model for a given 
value of the uncertainties and is given by = (^/ ) • 

3. Measurements 

3.J. Electrophysiological measurements 

The best way to access the electrophysiological quantities would be to directly mea¬ 
sure the electrophysiological potential at the heart surface. This is possible with an 
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Figure 3: Direct simulation results in long axis view — deformed mesh with interpolated transmembrane 
potential from electrical refined mesh. ^ 




















invasive procedure where a basket of electrodes is immersed in the ventricle and records 
the potential close to the endocardium surface m e.g.). In that case, the electrophysi- 
ological measurements can be represented by 

= (14) 


where the observation operator interpolates the values at the basket electrodes and 
X® gathers all the measurement errors. After spatial discretization 

Z®(t) =H®A® + x® (15) 


where now x® takes also into account the finite element discretization error. 

To avoid such invasive measurement procedures, it is preferable to rely on electrocar¬ 
diograms. Physiologically, the electrical potential diffuses from the heart to the rest of 
the body through the pericardium. The ECG depicts the time course of standardized 
potential differences on the body surface. We refer to m for examples of healthy and 
pathological ECG obtained with the electrophysiology model used in the present paper. 

The electrical diffusion within the body can be modelled by a Poisson equation 
with the extracellular potential as a boundary condition on the heart surface and homo¬ 
geneous Neumann boundary elsewhere: 

fy-(£^-V ut) =0, infl^ 

< mt = Ue, on SHq (16) 

[ (iZrp • y ut) • u = 0, on 


where the conductivity tensor is assumed to be isotropic, but non-homogeneous to ac¬ 
count for the different conductivities of the lungs, the bones and the rest of the body 
(see [in])- 

After spatial discretization - see the corresponding thorax mesh in Figure - the 
ECG measurements denoted by are still related to the electrical state variable by an 
observation operator i/® which is the composition of the discrete diffusion operator and 
a linear combination of the potential at 9 points of the body surface. Assuming that 
the model is accurate enough, the actual EGG can still be represented by (14) where 
X® gathers all the possible measurement errors, including the errors resulting from the 
modelling of the measurement procedure. After spatial discretization, a relation of the 
form (15) still holds, but with x® also accounting for the discretization error. 


3.2. Mechanical measurements 

The heart contraction can also be perceived from a kinematic and even a mechani¬ 
cal perspective. In fact, the displacements of the cavity can be observed using imaging 
modalities, for example Cine-MRI sequence, CT-sequence or even Tagged-MRI sequence. 
Gonsidering for example this last type of measurement and assuming adequate registra¬ 
tion strategies allowing to extract 3d displacements of material points [55j , we can assume 
that displacements are measured in a large part of the ventricles, typically between two 
short axis planes, the first one slightly lower than the base and the second one slightly 
higher than the apex - see Figure With more standard sequences, like Cine-MRI or 
CT sequences, we can benefit from optical flow strategies to extract displacements most 
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of the time on part of the boundary. In both cases, after discretization, we can assume 
the existence of an observation operator H“ such that the measured displacements are 
given by 




where y™ is now a vector taking into account the discretization error. 

Remark 2 (Generalized observation operators) 

Extracting displacements using registration techniques is still a challenge and it of¬ 
ten leads to the measurement of an apparent motion instead of the real motion. As 
a consequence, it can be more adequate to consider the structure measured in itself 
more than its collection of material points. This was done typically in |44) where in 
Cine-MRI and CT sequences the contours of interest are extracted and compared to 
the contours produced by the simulation by the use of a discrepancy field of the form 
= dist g^(x(£,t).t) = x(^,t). Hence, the closer to 0 the discrepancy is the more 
information we have extracted from the images. As demonstrated in |44j . this type of dis¬ 
crepancy operator is a simple non-linear generalization of the discrepancy computed by 
subtracting to the actual measurements the one produced by the use of the observation 
operator defined below. This generalization has demonstrated its efficiency in several real 
cases of data assimilation investigations mill]. But, for the sake of simplicity, and since 
this work is primarily aimed at providing a proof of concept, we limit ourselves to more 
simple observation operators. Note that the use of a non-linear observation operator is, 
in the end, always analyzed in the light of the corresponding linearized operator |44j . 

The stresses are more intricate to measure. For now, the most common measure¬ 
ment available is an intracavity pressure obtained after invasive catheterization or re¬ 
constructed from arterial blood pressure measurements. In this article we assume that 
we can consider the intracavity pressure as a given source term in our model Q. In a 
more general case, we could model the evolution of the intracavity pressures with the 
help for example of a Windkessel model and then use the possibly noisy measurements 
in a data assimilation context where the intracavity pressures are additional variables of 
the model. 

3.3. Multi-modalities measurements 

Finally the observations are concatenated and therefore given and the complete form 



(17) 


which for the sake of generality will not be necessarily considered as linear when it is not 
mandatory. 

4. Data assimilation principles 

Data assimilation has become a very popular strategy to estimate a wide range of 
modeling uncertainties in numerical simulations [5|. It has been initiated for environ¬ 
mental sciences but has now reached life sciences and especially cardiology. To introduce 
its main concepts in the state-space formalism, we consider a dynamical model 


X = A{x,d,t). 
12 


( 18 ) 


In this equation, x denotes the state variable, namely, the physical quantity which the 
model aims at describing during its time evolution. In this generic notation, the whole 
model is essentially summarized in the dynamical operator A, which applies on the state 
variable itself, and may depend on time t as well as on a set of physical parameters 
denoted by d. In this work, (18) corresponds to the electromechanical system dehned 
by (H and a, or - to avoid issues associated with PDEs - its discrete counterpart 
®) and ( |l^ (but most of what is presented here can be actually generalized to infinite 
dimensional systems HD- 


The initial condition a;(0) and the parameter vector d must be prescribed in (18). 
Hence, when these quantities are unknown, they have to be estimated. In general ^(O) 
and d are decomposed into known parts a;* and i?<> called a priori and uncertain parts 
and C’’: 

fx(0) = Xo 
\d = d^ + C^. 

When necessary, the trajectory of x will be denoted by 

Let us now re-consider an augmented state = (§) in order to integrate parameter 
uncertainties. The corresponding augmented model is defined by 


= 'dl('x), 


and ^a:(0) = 


Xo+C" 


= 'xo + 'C. 


(19) 


Hence the parameters and the state can be consider in the same formalism. However, 
these two quantities differ in their dimension. Indeed, when the state is defined based on 
a PDEs model, the dimension of the state variable after space discretization is typically 
of the size of 10^ to 10^ degrees of freedom. By contrast, the size of the parameter 
vector is generally much more limited. Even if distributed parameters are considered, 
their variation should be considered smooth enough so that they can be discretized on a 
coarse mesh or a subdivision of the domain into large regions. Therefore, typically less 
than a hundred of parameter values need to be estimated. 

The measurements presented in Section]^ can be cast in a general form: given a real 
trajectory x\ the noisy observations are represented using an observation operator H 
such that 

z = H{x\t) + x{l). ( 20 ) 

where z denotes the actual data field. Note that H can possibly be extended to a function 
of ^x and will then be denoted by 


4-.1. Data assimilation by filtering 

A common strategy in data assimilation is to minimize a criterion that balances the 
confidence in the a priori value of the state and parameters and a discrepancy measure 
between the given observations and the simulated ones. In the state-space formalism this 
criterion is typically least-square and reads 

^(C) = ^ll'CII^-1 + I l^\\z - ^Hi\.^^)\\l,dt. (21) 

where M is a metric on the observation space and Po is the inverse of a metric which 
can be interpreted as an initial uncertainty covariance. The so-called 4D- Var method 
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|33| consists in finding the unknown quantities '’C by minimizing ^ under the constraint 
of following the dynamics ( |18[ ). When T goes to infinity the minimization is expected 
to produce a better and better estimate. This minimization under constraint requires 
the computation of an adjoint variable used to compute the gradient. Hence, a gradient 
based descent algorithm requires numerous iterations of the direct and adjoint dynamics. 

When we do not wish to precisely retrieve the initial error but only seek to ac¬ 
curately approximate the state “independently” of the possible initial error, 

we can avoid minimization iterations by the use of sequential estimation methods. The 
principle in sequential estimation methods is to introduce a modified system - denoted 
by a circumflex accent - called observer whose dynamics is changed to incorporate a cor¬ 
rection based on the measured discrepancy. The new system dynamics therefore reads 

= A{^x, t) + G{z - '’H(’x)), ^2;(0) = '’x^ (22) 

where G is called the observer gain or filter. The ultimate objective of the observer ^x{t) 
is to converge in time to the real trajectory 


and the gain has to be designed regarding this objective. Two strategies are commonly 
followed in this respect. First, the gain can be constructed from the optimality criterion 
(21) by defining the so-called optimal ofeserwe?!^- or optimal sequential estimator - with 

^[argmin 


An optimal gain is obtained by differentiating this definition with respect to the time 
variable t appearing in both the trajectory ^x{t) and in the criterion This ap¬ 

proach is well-known in a linear framework - i.e. where all the operators are linear 
- and leads to the famous Kalman-Bucy filter [11201 EH- The optimal filter is then 
defined from the solution of a Riccati equation. In a non-linear framework, the gain is 
more intricate to compute and derives from a Hamilton-Jacobi-Bellman solution |23j . 
Hence, in this case, numerous works rely on approximate solutions as defined for exam¬ 
ple by the Extended Kalman Filter (EKF) which uses the Riccati equation of the linear 
Kalman filter with the tangent operator of the non-linear model and observation oper¬ 
ators m- Eventually, the great advantage of the optimal sequential approach is that, 
like the 4D-Var, it can be defined for every model and every observation operator. The 
biggest drawback is that - even with approximate solutions - the filter is very costly to 
compute, especially for large dimensional systems like the ones produced by PDEs or 
their discretizations. Therefore, alternative strategies - the Luenberger observers and 
the Reduced Order Optimal filtering - will be presented below. 


4-2. The Extended Kalman Filter (EKF) 

Even if we will have to rely on alternative strategies, it is interesting to develop the 
optimal filtering equations. For the sake of simplicity, we proceed after space discretiza¬ 
tion. The space discretized system variable is denoted by X and the space discretized 


1 


the observer is called optimal due to the fact that it is associated with an optimal control problem. 
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parameters by 6. The observer state X and parameters 9 follow a modified version of 
the dynamics: 

= A(X, 9) + G^{Z- H(l,t)), 

/9 = G%Z-Y{{X,t)), (23^ 

A(0) = Ao, 

9{Q)=9^, 

where the gain G = (^ ) is a linear operator. In the augmented form with ^X = (A) 
we write 

= ^A(X) + G[Z - ^(X, t)). 

The most classical gain is given by the Extended Kalman Filter (EKF): G = P(dH)'''M, 
where P is obtained from the solution of a Riccati equation on the augmented form 

P = (d^A)P + P (d^A)T - P(d^)TM(d‘H)P. (24) 

The operator P solution of the Riccati equation is called covariance since it can be linked 
in a stochastic framework to the covariance of the state estimation error evolving during 
the sequential estimation m- When decomposed on the state and parameter 

pxx px0\ 

this gives 

Qx ^ G® = (P^^)T(dxH)TM, 

and 

pee ^ _(pxe)T(d^H)TM(d^H)P^®, P®®(0) = P, 

p^» = (dxA)P"^® + (deA)P®® - P"^^(dxH)TM(dxH)P-^®, P^«(0) = 0 

p^^ = (dxA)P'^^ + (deA)(P^®)T + P^^(dxA)T + P-^®(deA)T 

-P"^^(dxH)TM(dxH)P"^^ P^"^(0) = Po 

(25) 

The practical algorithms - for instance time-discrete EKF or UKF for Unscented 
Kalman Filter - derived from this formulation are presented in [Appendix A[ Even 
if UKF differs from EKF, it is easy to prove that their analysis which relies on the 
linearization of their respective estimation error is based on the same error system |42j . 
This is due to the fact that when all operators are assumed to be linear, a finite difference 
operator or the tangent operator are in fact identical. This implies that we will rely on 
UKF for the numerical simulation whereas, for the estimator analysis, we will keep EKF 
which can easily be written in both time-continuous and time-discrete formulations. 

Various sequential estimators can be defined on any dynamical system, and parameter 
identification can be easily included. This is due to the underlying optimal principles 
behind the estimator presented before, even if in a non-linear framework we have only 
relied on approximation rules to extend the optimal filters found in the linear case. 
However, these estimators have a fundamental drawback since the covariance operator P 
is a full matrix of MjVx+d(lEf)- Recalling that is the dimension of the state variables 
which discretize the PDFs field variables, the optimal estimators are intractable for 
classical finite element models. 
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4-3. Reduced-Order Optimal Filters 

A classical first strategy for circumventing the curse of dimensionality consists in 
assuming a specific reduced-order form for the covariance operators. For example, making 
the ansatz 

Vt, P(t) = L{t)Vit)-^h{ty (26) 

with U an invertible matrix of small size r and L an extension operator, we can show |4dj 
that with linear operators the solution of the Riccati equation (24) in augmented form 
reduces to 

L = ^AL and U = LTRTMHL, (27) 

which are now computable in practice. In a non-linear framework, the covariance dy¬ 
namics can then be approximated as in |4dj by extending (27) as 


L = (d^A)L and U = (dH)TM(dH)L. 


(28) 


This observer is called Reduced-Order Extended Kalman Filter (RoEKF). The time dis¬ 
crete version of this approach is presented in|Appendix B[ as well as its UKF counterpart. 


4-4- Luenberger Filters 

Another way to circumvent the curse of dimensionality is to built a filter which is 
not based on an underlying optimal criterion. This idea was initially introduced in |35j 
and therefore is often called Luenberger filter and observer. It was quickly popularized 
in data assimilation [261I5811I] where the curse of dimensionality was very limiting for 
systems coming from the discretization of PDFs. In data assimilation this strategy is 
also referred to as the nudging approach because the filter is designed to “gently” correct 
the dynamics. The goal of the correction is to make the estimation error 
converge to 0. In a linear framework, the state error dynamics is 

= (A - '’GHfi, 

thus the strategy is to design a filter gain G based on the underlying physics of the 
problem which makes the dynamics operator {A — GF[) dissipative - strongly if possible. 
Such Luenberger filters are designed for the state estimation of specihc physical systems 
and for specific observation operators. For instance, filters have been proposed for trans¬ 
port equations Schrodinger problems [5], for waves m, beams, elasticity [IHIHI] or 
fluid-structure problems [5]. In Section [5.2[ this approach will be used for the mechanical 
problem. 


4-5. Joint state and parameter estimation and coupled system estimation with Luenberger 
filters 

System (19) has been defined from the augmented state vector gathering the original 
state and the parameters. After decomposition on x and i?, this implies that the initial 
parameters now evolve in time within the observer 


X = A{x,d,t)-\-G^{z — H{x)), X = Xo, 
h = G^{z-H{x)), i?(0) = do, 


( 29 ) 
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We have seen that the optimal filter does not make any difference between state and 
parameters. However, when dealing with Luenberger observer, there is no direct way to 
incorporate the parameter identification - i.e. to design - from an already defined 
Luengerger filter G® on the state. This problem is known as adaptative filtering O Ibbj . 
In particular in [42l [43] a systematic strategy has been proposed to extend a possible 
Luengerger observer to a joint state and parameters estimation by combining the Luen¬ 
berger filter and an optimal filter reduced to the remaining parameter space. Taking into 
account that the parameter space is of much smaller dimension, this strategy allows us 
to compute a physically-based gain on the large dimensional state and an optimal-based 
filter on the small parameter space. The resulting observer is proved to converge under 
adequate assumptions in linear cases and can be extended to non-linear cases. 

When considering multi-physics coupled systems, as the electromechanical problem 
(HI), the same type of problem occurs as for joint state and parameters estimation. If 
optimal filters are considered, there is no difficulty to combine them. When measurements 
are available, and when Luenberger filters have also been designed for the different types 
of physics, it is also easy to combined them. But, as said before, optimal filters are 
extremely expensive, and it is not always possible to design a Luenberger filter for all the 
physical compartments of the problem. It will be justified in Section |5.3| that, as done 
for the state and parameter estimation, it is possible for one-way coupled problems to 
combine a Luenberger filter in one part of the system and a reduced optimal filter on the 
other part. 

5. Data assimilation for the electromechanical model 

5.1. State estimation in electrophysiology: a reduced-order approaeh 

The above data assimilation principles can been applied to the system modeling the 
cardiac electrophysiology. Using a few measures of the electrical potential, the goal is 
to reduce the uncertainties on the state and the parameters involved in the modeling of 
cardiac cells. For this system, contrary to the mechanical system that will be presented in 
the next section, there is no straightforward Luenberger observer. The state has therefore 
to be filtered differently. 

A possible strategy to overcome the curse of dimensionality induced by the optimal 
filtering is to discretize the problem on a low-dimensional basis. Here, we propose to build 
this basis by Proper Orthogonal Decomposition (POD). The POD basis is obtained by 
keeping the most relevant modes resulting from a Principal Component Analysis of a set 
of pre-computed solutions. For more details about POD, we refer for example to |311|S3] 
and to [m [16] for applications to electrophysiology. 

We denote by H € Mjv'.atpqd(®-) the matrix made of the first A^pod modes. These 
modes are orthonormal with respect to a given scalar product, typically in 
or Denoting by Mpoo the Gramian matrix associated with this scalar product, 

the POD expansion coefficients of a vector X® are given by a = H^MpodA'®. Defining 
XpoD = Ila and XJ_ = X® —Ha, the state vector can be decomposed as X® = X®Qp,-|-X^. 

As for the reduced filtering, a first strategy would consist in replacing the full order 
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model by the reduced one: 


{ ve 
^POD 

^POD 
0° = 


(0) = + Cx', 


+ Cs” 


where Ap^p is computed from the finite element matrices M®, Kf and the vectors F®, G 
are projected on the POD space. Then, after discretization, this (small) system can be 
estimated by EKF or UKF (Appendix A|. In doing so, the estimator is known to be 
stable, but this approach has a drawback: if the POD basis is not rich enough to capture 
the relevant phenomena, the reduced order dynamics may be a poor approximation of 
the full order one. Thus we may have a consistency problem in our way of modeling 
and approximating the system of interest. That is why we prefer another strategy which 
consists in solving the full order dynamics, and then applying the filter to the projection 
of the state vector on the POD basis. In other words, we apply RoUKF algorithm 
(Appendix B) with a reduced variable made of the parameters and the components of 
the augmented state vector on the POD basis: 





Jf^iVpOD+d ^ 



(30) 


Contrary to the first approach, the full model is preserved. The drawback is that it 
cannot be proved that the part of the error which is not filtered will not grow and then 
pollute the approximation. Thus, while the first approach could suffer from a consistency 
problem, the second one can suffer from a stability problem. Nevertheless, in all our test 
cases, this second strategy proved to be robust and accurate. 

By combining this strategy with UKF, the following algorithm is obtained to filter 
the electrical state. The initial condition projector is given by 


L(0) = 





Then we consider an adequate UKF sampling rule composed of weights and particles (see 
[29] or Appendix A), we store the associated weights (ai) in the diagonal matrix Da 


and precompute specific particles, here the so-called unitary sigma-points (i.e. with zero 
mean and unit covariance). More specifically, we consider the p = iVpoo + d + 1 unitary 
simplex sigma points, where iVpoD -I- d is the reduced space dimension [551135] . We denote 
them by (/*)i<i<p and perform at each time step: a first sampling step which generates 
the particles identified by a subscript a prediction step denoted by an additional 
superscript “ and a correction step denoted instead by an additional superscript 

1. Sampling: 


Cn 



^±n 

= + lV • • /* 

, l<i< 


= d++D„•C^/^ 

1 < i < p 

t'n 

= e+ + Li-ci-p, 

1 < i < p 


(31a) 
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2. Prediction: 


3. Correction: 



/■ ^ 

w- 

n+1 

— A„_|_i|„(x|t^~'’), 

1 < i < p 



.[i]- 

<+1 

= nTl[:;-, 1< 

i < p 


< 

1 

N- 

±n+l 

= x'Vi - nd'ti^ 

1 <i < P 




±n+l 

— tW- 




b^n+l 

= I.i=i aia«+i 





+1 

= ELi«4ti 



' T^_L 
^n+1 

= 

\.^±n 




I^n+1 

= 

(“n+l 




-^n+1 

= 

V^n+l 

]^?a[/bl]T 



^n+1 

= 

Hn+1 




^n+l 

= 

y-P 




Tn+l 

= 

fyl*]- 




U«+1 

= 

i + r);+iW-|ir„+i 




= 

^-Ln+l “b bn+lMn+iFj^^g 

M„+i(Z„+i 

~ ^n+l) 

b^n+1 

= 

b^n+1 


n+1 (-^n+l 

^n+l) 


= 

^n+1 


n+1 (-^n+l 

^n+l) 


(31b) 


(31c) 


where [/bl] is the matrix concatenating the (/*) vectors side by side, and similarly for 
other vectors [H]. From the last three corrections defined in (31), the state correction 
reads 

1++1 = 1-+1 + + nE;+0 u"+ir^+iM„+i(z„+i - z-^,). 

5.2. State estimation in mechanics: Luenberger observers 

To filter the mechanical state, a Luenberger approach is used, following [44] . Consider 
the first two equations of © with 


dty = V, in 

pdtv — V • (T) = 0, in fig 


(32) 


where the observations are given by 


^ =2/i, 


(33) 


We then define the estimator by 


'dty = v + 7Ext„^b.(z“ - 


J, 


pdtv-Y-{T)=0, 

^same boundary conditions as ([^ 
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in fig 
in fin 


(34) 





where 7 is a scalar gain and Ext^obs is an extension operator typically given for any 
displacement field d by 


^cxt 


Ext^obs (d) 


'V-(4:£(d®"*)) = 0 
= d 

A - ■ n = ksd^"\ 

{A : £{d±'^^)) • H = 0, 


in 

. ,obs 

m LO^ 

on r„ 

on 


(35) 


with A the elasticity tensor coming from the linearization of T around 0 or a given 
trajectory at a given time t. On a linearized system, it can be proved that the the state 
error x = (|) = ( 7 _|) between the estimator and the target tends to zero |44j . 

A consistent space discretization is given by 


Y = V + 7Ext(Z“ - H“(A)), 
M“y+ K“(y,E) = N“ 


(36) 


which converges to the real trajectory. Note that the proof of convergence after dis¬ 
cretization is in general difficult and should be proved here using the internal viscosity 
to control the spurious high frequencies introduced by the discretization. 

Finally in a state-space form, there exists a Luenberger filter G™ such that 

i“=A“(l“)-hGj(^(Z”-H“(A“)) (37) 

which converges for any initial error to an observed trajectory. This Luenberger filter, 
defined in the context of passive non-linear mechanics, was also shown to be robust to 
the introduction of the active part in the cardiac mechanics mill]. In this context, the 
Luenberger filter can be considered as a filter on the reduced space of the displacement 
and velocity field whereas the internal variables associated with the heart contraction 
are already stable. Nevertheless, note that the filter robustness has been demonstrated 
numerically but remains to be proved theoretically. 


5.3. Estimation of one-way coupled systems 


In this section we present an original strategy to aggregate the already defined filters of 
each submodel - namely the electrophysiological model and the mechanical model. This 
strategy is then justified by an estimation error analysis. To simplify the presentation 
and analysis of our method, we first consider the following one way coupling system 
with perfectly known parameters. Indeed, we have already seen how a joint state and 
parameter estimation can be added “in a second stage” once the state estimation is 
proved to be effective. Hence in practice the parameter identification will be considered 
but here, we can focus on the state as the main difficult aspect of the analysis. We thus 
consider 


A® = A® (A®), A®(0) = A® -h C 

A*" = A“(A“,A®), A“(0)=A“-hC“ 


(38) 
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and we propose to combine the Luenberger state filter on the mechanical part and the 
reduced order filter of Section [STT] on the electrical part. We also recall that the observa¬ 
tion operator is composed of two concatenated operators H®(X®) and H“(X“) as defined 


in (17). 


We can consider a RoEKF without loss by generality at the time-continuous level. 
Concerning the derivatives, we denote dm and de respectively the partial differential with 
respect to X™ and X® respectively, whereas we keep d when there is no ambiguity or 
when the differential is total - i.e. with respect to all variables. From the general reduced 


order formulation of the FKF presented in Section 4.3 we introduce this time 


L = 


such that 

'X® = A®(X®)-tL®U-iLTdHTM(Z-H(X)), X®(0) = X® 

X“ = A“(X“,X®)-tG£'(Z“-H“(X“)) 

-bL“U-iLTdHTM(Z-H(X)), X“(0)=X“ 

< L® = (dA®(X®))L®, L®(0) = 1 (39) 

L“ = (dmA'"(X“, X“) - G5"dH“(X'"))L“ 

-b(deA“(X“, X®))L® L'"(0) = 0 
V = LT(dH)TM(dH)L, U(0)=U|. 


Note that in the dynamics (39) the Luenberger filter only applies on the mechanical part 
using only the mechanical data. However, the optimal filter strategy allows to benefit 
from both the electrical and mechanical data to correct the electrophysiology dynamics. 
This correction is then reverberated to the mechanical dynamics which also reads 


X“ = A“(X“, X®) -b Gf^{Z^ - H“(X”)) -b L“(L®)-i(X® 


A®(X®)), 


when assuming L® invertible. Note that in the particular case where A® is linear the last 
expression simplifies into 


X“ = A“(X“, X®) -b G5^(Z“ - H“(X“)) -b L“ (L®)-iX®, 

5.3.1. Convergence analysis 

To ensure the convergence of the observer, L® is assumed to be invertible. Note 
that L® follows the tangent dynamics of the electrophysiological around the estimated 
trajectory of the model, starting from the initial condition L®(0) = 1. Therefore, even for 
a dissipative system, we can assume L® to be invertible but potentially ill-conditioned. 
Moreover the invertibility of L® can also be verified numerically in the specific case of 
concern. We recall that we denote by X = X^ — X the estimation error. The variable X 
follows a non-linear dynamics. Then by linearization of this dynamics around the target 
trajectory we get a linear dynamics satisfied by linearized estimation error denoted by 
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( 40 ) 


5X = ( ). In fact we have 

' dX ' 


{ 5X<^ =-U-iLT(dH)TM(dH)«, J1«(0) = C® 

= (d„iA - Gf dmH“)<5A“ + (doA)<£¥®+ 

L“(L®)-MA®, dA“(0) = C“- 

Here we introduce the change of variables 

(<5A“, 6X^) (Jr?, 5^l) = (<5A“ - L“(L®)-MA®, (L®)-MA®), 

and obtain 

r Sr, = (d^A - G^d^H“)5ry, Sr,i0) = ^ 

}s,i =-U-iLT(dH)TM(dH)L(5Ai 
[ -U-iLT(dmH“)TM“(dmH'")(5r7, (5A®(0) = C- 


(41) 


The first equation corresponds to the dynamics of the linearized error studied for the 
mechanical system. Hence it converges to 0. Therefore the second term in the second 
equation tends to 0. The homogeneous part of the second equation can then be proved 
to converge to 0 if the following observability (42) condition is satisfied with our linear 
observation operator ~ namely, dH = H in our particular example. Namely we expect 
for all initial error dA®(0) that 


3(C,T), r\\R{L{L^)-^SX^)\\l,>C\\dX^m%, (42) 

Jo 

which can be at least verified numerically. In the last observability condition we see that 
L(L®)“^dA® represents the effect of a variation on (5.A® on both the electrophysiology 
and the mechanics. This effect is then observed through H. The observability is thus 
expected to be improved with respect to the situation where only the electrophysiology 
is considered. This will be confirmed numerically. 


5.3.2. Practical algorithm 

The correction steps of algorithms (B.5) (for mechanics) and (31) (for electrophysi¬ 
ology) are independent of the problem, and can thus be implemented by adopting again 
a master-slave strategy. In our electromechanical simulator, this yields a two-level in- 
terweaved master-slave strategy, the master handling the electromechanical problem be¬ 
coming a slave of the sequential estimator algorithm. The tasks of the electromechanical 
master are the data exchange and the correction of unknowns and parameters. The 
estimation algorithm is summarized in Figure]^ for a generic time step. The first data 
transmission appearing in Figure is done in order to update the initial value for each 
solver. This is handled by the electromechanical master. 


6. Numerical illustrations 

In this section, the methodology presented in Section [5T] is applied to the electrophys- 
iological problem of Section |2.1[ The purpose is twofold: first, assess the reduced order 
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Figure 4: Estimation algorithm: “particle solver” denotes the electromechanical solver presented in 
Section EH 


state filtering based on POD; second, illustrate the interest of estimating simultaneaously 
the state and the parameters. 

All the test cases of this Section are performed on the geometry of Figure including 
the fibers in the conductivity tensors, and with the parameters typically used to generate 
healthy EGG [TU]. Synthetic data are generated by applying in Q an external current 
lapp during 25 ms. Two sets of measurements are considered: (1) 12-lead EGG only or 
(2) 12-lead EGG enriched with the extra-cellular potential at 8 epicardium points. The 
first set of measurements is of course the easiest to obtain. The rationale for the second 
set is to foresee what could be the benefit of including the measurements available from 
cardiac stimulators or implantable defibrillators. For both cases, these measurements 
are perturbed by an additive Gaussian white noise of standard deviation of 0.25mF 
(Figure [^. Note that another option could be to include endocardial extracellular po¬ 
tentials, if a catheter can provide them. We do not anticipate significant differences with 
respect to the results obtained with epicardial measurements. 


6.1. ECG based state estimation of the electrophysiological model 

Before addressing parameter estimations, a first test is run to assess the reduced order 
state filter, assuming that all parameters are perfectly known - for instance with Tin = 0.8 
and Tout = 18. In other words, we study the problem of the state estimation alone applied 
only tp the electrophysiological model. The POD Basis used in the reduced order state 
filter (31) was obtained, in this case, from snapshots generated with these parameters, 
among others. We choose to run the estimation simulation without any external current 
and compare it with the target simulation starting from t = 40 ms. When the estimation 
starts, the state variables are therefore significantly perturbed with respect to the direct 
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simulation used to generate the synthetic data. In doing so, we introduce a significant 
state error in our estimation problem. Figure compares the estimated mean value of 
the transmembrane potential Vm and the ionic variable w with the target values, while 
Figure compares the spatial distributions over the whole heart domain. We see that, 
despite the very large initial error, the estimator succeeds quite well in compensating for 
the lack of information about the initial stimulation and retrieving in time the accurate 
electrical state. Indeed, it is remarkable that with a correction only coming from the 
reduced state filter and applied from t = 40ms, the resulting estimated transmembrane 
potential and ionic current are in good agreement with their target values. Eventually, 
we point out that, as expected, adding a few measurements on the myocardium slightly 
improves the result. However this type of measurement clearly requires an invasive 
procedure m- 


Simulation 

1 

2 

3 

4 

En 

0.8 

0.8 

1.2 

1.2 

'^out 

18.0 

12.0 

18.0 

12.0 


Table 1: Set of parameters used to build the POD basis. 



-Direct Simulation 

-EGG based Estimation 

-EGG + 8 epi. pts based Estimation 


Figure 5: Estimation of the transmembrane potential Vm (Left) and the ionic variable w (Right). The 
direct simulation was performed with an external current applied during 25ms. The curves represent 
space-averaged quantities. 


6.2. ECG based parameter identification of the electrophysiological model 

Next, we address the question of the identification of parameters Tin d-nd Tout when 
still only considering the electrophysiological model. The values used to generate the 
synthetic data are Tin = 1 ^tid Tout = Id. The initial guess or a priori in the inverse 

24 












Figure 6: Space distribution of the transmembrane potential for the direct simulation (top) and the state 
estimation with ECG measurements (bottom). There is no initial activation in the bottom simulation 
in order to generate a strong error in the initial condition. The correction of the filter is applied from 
i = 40 ms. Note the good agreement of both simulation from t = 55ms. 


Tin estimation Tout estimation 



Parameter Identification only 
Parameter Identification Standard Deviation 
-Target Value 

-Joint State-Parameter Estimation 

Joint State-Parameter Standard Deviation 


Figure 7: Comparison on rin (Left) and Tout (Right) estimation, with (black line) and without (gray line) 
reduced state filtering. These curves clearly show that reducing the uncertainties on the state actually 
improves the identification of parameters 
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problems is Tin = 1-5, Tout = 11- We point out that neither the values used to generate 
the synthetic data nor the initial ones are included in the set of solutions used to generate 
the POD basis. In fact, the POD basis used in the reduced filter (31) is computed from 
snapshots obtained with 4 sets of parameters Tin and Tout “ presented in table - which 
are in a neighbourhood of the expected value, but do not coincide with it as theoretically 
justified in |16j . Before launching the estimation, we proceed to a reparametrization of 
the form Tin = 2 ^^T^n where and T^nt are given and 6 i and 62 

are the new values to be estimated. This is motivated by the fact that the parameters 
should be maintained positive during the estimation and the uncertainty variance is more 
naturally centered with respect to a power of 2 of initial parameter. We then compare 
two strategies: (1) parameter estimation only; (2) joint state-parameter estimation. In 
strategy (1), the component to be filtered is limited to the parameters 9, whereas in 
strategy (2), X'^ includes the parameters 9 and the expansion coefficient a on the POD 
basis (see (30)). 


Param. 

Target 

Value 

A priori 
(%Error) 

Only Param. 
Identification 

Joint State & 
Parameter Estim. 

Tin 

1.0 

1.5(50%) 

2.19(118.8%) 

1.02(1.99%) 

Taut 

16.0 

11.0(31.25%) 

6.96(56.48%) 

15.15(5.3%) 


Table 2: Target vs A priori values for rin and Tout 


The evolution of Tin and Tout is reported in Figure and summarized in Table In 
particular we clearly see in Figure the convergence of the estimator on the parameters 
variable. As an additional illustration of the performance of the joint state-parameter 
estimator, we plot in Figure [^8 leads of the corresponding standard body surface ECGs 
obtained with the initial guess of the parameters (curve named “EGG built from wrong 
initial guess”), with the parameters obtained after a simple parameter estimation (curve 
named “EGG after identification only”), with the estimated parameters obtained after 
a joint state-parameter estimation (curve named “EGG after joint state and param. 
estim.”). All these curves should be compared to the one named “Target EGG” which 
corresponds to the measurements used for the estimation. Two comments are in order. 
First, we note that the EGG obtained with the “wrong initial guess” is indeed very 
different from the target. This means that our initial guess was far from the solution we 
looked for. Second, we see that the joint state-parameter estimation clearly outperforms 
the simple parameter estimation (as noted in a different context in [43]). This means that 
only a joint state-parameter estimation allows us to produce an electrical state that is 
really compatible with the observations starting from realistic initial covariance - namely 
without forcing a priori the estimation to fit the data. 

6.3. Full electromechanical data assimilation 

In this section, we demonstrate numerically the efficiency of our complete state and 
parameter estimation chain for coupled models and, more importantly, we show that the 
estimation of a piecewise constant electrical parameter can be improved by enriching the 
electrical observations with kinematical observations. The parameter of interest is here 
Tciose, which controls the plateau duration of the action potential of the cardiac cells and 
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Figure 8: ECGs after state and/or parameter estimation (zoom on a time window of 400 ms) 
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is supposed to be heterogeneous in the myocardium. In fact, Tdose is assumed to take 
different constant values in 4 regions: inner part of the myocardium in the left ventricle, 
that will be called “endocardium” for simplicity; outer part of the myocardium in the 
left ventricle, that will be called “epicardium”; a thin region between the endocardium 
and the epicardium in the left ventricle, that will be called “M-cell”; and finally the right 
ventricle. 

Two kinds of synthetic observations are generated from a direct simulation: the 
electrocardiograms and the displacements of the myocardium. Hence for this numerical 
illustration, we have at our disposal a 3d field of displacement as it could ideally be 
processed from a 3d tagged-MRI sequence [55]. However this is clearly an ideal situation 
from the mechanical point of view. Nevertheless, we believe that this illustration offers 
a clear insight into the maximum of information that could be earned with mechanical 
observations from the point of view of the electrophysiology model. 
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Figure 9: I lead electrocardiogram with additive Gaussian noise (left) and displacement of a mesh point 
with additive Gaussian noise (right) 


These quantities are perturbed by an additive 10% gaussian white noise in time and 
space, proportional to the mean amplitude of the signal - see the corresponding signals in 
Figure]^ This type of noise is clearly illustrative as it is reasonable to choose independent 
time observations whereas the spatial distribution of the noise could certainly be more 
complex than a gaussian noise. In fact, real observations can even be spatially biased. 
However, this choice of noise is well adapted to evaluate a minima the robustness to noise 
of an estimation method. Facing these noises we follow the recommendation made in 
m with real data to choose a time-discretized observation norm built from the inverse 
of the standard deviation Jobs by 

Vn,M^ = Ate(J:bs)-"l, and ,M“ = At„(J“ 


with ~ 10% * 2.5 mV = 0.25mV and ~ 10% * 10 mm = 1mm. 


■^obs 

As in Section 


6.2 


the unknowns are reparametrized as Tdose = 2 




close’ 


where 


close 


corresponds to the initial guess, 
gathered in Table 


The values to be estimated and the initial guess are 
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Figure 10: Tdose estimation on endocardium, m-cell, epicardium and right-ventricle regions, comparison 
ECG observations vs ECG+mechanical observations 
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Region 

Target value 

Initial guess 

Error(%) 

Endocardium 

140.0 

56.0 

60% 

M-cell 

105.0 

42.0 

60% 

Epicardium 

105.0 

42.0 

60% 

Right ventricle 

120.0 

48.0 

60% 


Table 3: Target value and initial guess of Tdose in ths different regions 


ECG I Lead Sensitivity 



f 


Mech. Sensitivity 
■ 10 -'* 


0.5 


0 

-0.5 

-1 


RV 

^^ ' 

200 400 600 800 


Time (ms) 



Figure 11: Measurements sensitivity with respect to Tdose- ECG Lead I sensitivity (left) and mechanical 
sensitivity (right, the four curves are superimposed) 


Figure]^ shows the time evolution of Tdose during the sequential estimation in each 
region for electrical measurements only and for the electromechanical - ECG plus dis¬ 
placements - available measurements. The benefit of the electromechanical measure¬ 
ments is striking, especially in the M-cell region. Then we present the hnal estimated 
values and the error with respect to the target values in Table We also report the 
final estimated values for electromechanical measurements after three heart beats where 
the estimation error is even more reduced. With the electrical measurements only, we 
were not able to run three heart beats, because the estimated parameters diverged too 
much from physiological values. On the contrary, with electromechanical measurements, 
the results keep improving along the three heart beats. Thus, it clearly appears that the 
estimation is much more accurate and robust when the electromechanical measurements 
are taken into account. Our data assimilation procedure results can also be evaluated in 
the light of ECGs that the model can produce. Indeed, we plot in Figurej^the reference 
ECG (in dashed black), the EGG corresponding to the initial guess (in black), the EGG 
obtained from the estimated values with the electrical measurement (in gray) and finally 
the EGG obtained from the estimated values with the electromechanical measurements 
(in dashed gray). We clearly see that the EGG corresponding to the parameter initial 
guess has a completely false T-wave, whereas after the estimation procedure the ECGs 
are much closer to the reference, especially with electromechanical measurements. 

In order to quantitatively understand the relevance of using the mechanical observa- 
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Parameter 

EGG only Estim 

EGG-|-Mech Estim 

EGG-|-Mech Estim 


after 1 beat 

after 1 beat 

after 3 beat 


(Error %) 

(Error %) 

(Error %) 

,-endo 

' close 

121.63 (13.12%) 

139.42 (0.41%) 

140.9 (0.64%) 

—niccll 
'close 

70.12 (33.22%) 

96.98 (7.64%) 

104.23 (0.73%) 

' close 

94.64 (9.87%) 

102.68 (2.21%) 

104.28 (0.69%) 

^RV 

close 

106.94 (10.88%) 

116.26 (3.12%) 

118.74 (1.05%) 


Table 4: Identification of Tdose in ths various regions and for the different scenarios 


tions to estimate an electrical parameter, a sensitivity analysis is performed. To simplify 
the presentation of this last calculus, let us consider a time continuous ROEKF for the 
parameter estimation. We recall that the sensitivity of the model with respect to its 
parameters is directly estimated with U in the ROEKF filter |33]. In our practical case 
we recall that IT is a matrix which consists of 4 columns vectors (one for each region 
of the Table of dimension the number of the electrical model degrees of freedom plus 
the mechanical model degrees of freedom. Therefore, the measurements sensitivity with 
respect to the parameters is estimated by HIT. Note that we have already seen such 


sensitivity terms in the observability condition (42) introduced in the convergence study 
of our coupled estimator. In our practical case, we focus on the sensitivity s®, resp. s™, 
of the first lead of the ECG, resp. the mean value of the measured displacements norm 
in the myocardium, with respect to Tdose, and we consider normalized quantities. We 
denote by {Z’^)i the first index of Z® where the hrst lead of the ECG is gathered. We 
thus define ^ 

s%t) = !T>^(H®L®(f))i, with = U \{Z^{t)U dt, 

Zi 4 Jo 

and T is typically the heart beat duration. Identically we dehne 


i/,\ Tcloseo 

' ' zjm. 


((H“L“(f))TM“H“L“(<))5, with z’^ = ^ f dt. 


In practice when using the ROUKF the computations are very similar but time in¬ 
tegrations is replaced by time iteration summations whereas HL(t) are replaced by 
computed from (31). Our two normalized sensitivities are represented 


in Figure [m for the four different regions. In the epicardium and in the right ventricle, 
the sensitivity of the electrical measurement is higher than the mechanical one, as ex¬ 
pected. But interestingly, the electrical sensitivity is lower than the mechanical one in 
the M-cell and endocardium regions. This is particularly the case here where we observe 
the displacement in the whole heart. However recent advances in tagged-MRI allows to 
expect that intra-myocardial kinematics will be available in the future whereas it is out of 
reach for non-invasive electrical measurement. In addition, the mechanical observation 
is affected by Tdose over a longer time window. Hence, the mechanical measurements 
therefore increase the time interval during which the observer has the opportunity to 
correct the value of Tdose- 

In conclusion, this last example is a rich illustration of our sequential strategy, with 
a state filter based on a Luenberger approach for the mechanics and on a reduced order 
optimal filter for the electrophysiology and the parameters. In particular it allows us to 
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Figure 12: ECGs obtained before and after estimation, compared to the ECG used as a measurement 
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numerically justify the importance of having a correction for every type of uncertainties, 
especially for state uncertainties. Moreover we numerically demonstrate the interest of 
taking advantage of the multi-physics nature of a problem to identify some parameters 
“ here Tdose accounting for the plateau duration - which have an impact on the coupled 
physics. 

7 . Conclusion 

In this work, a complete framework for the joint state and parameter estimation of a 
complex multi-physics problem has been presented. The strategy consists of Luenberger 
observers and reduced order optimal observers. The robustness of the approach is ensured 
by the fact that each component of this system is corrected by the data. In particular, 
the importance of state estimation when performing a parameter identihcation has been 
illustrated. For weakly coupled systems we have justified that our hlter combination and 
aggregation lead to a reduction of the estimation error. It is expected that the same 
type of results could be proved for fully coupled problems in the future. This general 
strategy has been applied to an electromechanical model of the heart but other coupled 
systems of interest in various fields of application can be estimated with this approach. 
In this framework, we have also shown that the use of mechanical data could considerably 
improve the observability of the problem. 

The direct electromechanical model is clearly more demanding than the direct elec¬ 
trophysiology model. However, we have shown that the associated inverse problem which 
consists in identifying its parameters from all the available data - mechanical and elec¬ 
trical - is better defined than the pure electrophysiology inverse problem. 

In conclusion, this work provides a promising new strategy to address the classical 
inverse problem of electrocardiography, and suggests a new way to reduce its well-known 
ill-posedness. The numerical results, based on synthetic data, are presented primarily 
to validate the methodology. This validation being successfully achieved, an important 
issue to be addressed in future works concerns the adequacy and accuracy of the models. 
Another important perspective is the validation against real data. Using real data will 
necessitate a complete patient workflow where anatomical data, EGG and Tagged-MRI 
(or at least Gine-MRI based on the results of M) are acquired and post-processed. In 
this respect, we point out that our strategy is indeed quite demanding in term of data 
acquisition and processing. This must be ultimately compared - in term of estimation ro¬ 
bustness and accuracy versus difficulty- to existing methods such as electrocardiographic 
imaging (EGGI) [H] where only EGG-like measurements are necessary with special vests 
recording up to 224 body-surface potentials. 
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Appendix A. Time-discrete EKF and UKF filters 

The time-discretization of the optimal estimator presented in Section |4.2| is based on 
the principle that the optimality should be conserved at the discrete level. Let us first 
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denote a stable and consistent discretization of the original model by 


^n+l — ^n) 

^0 = + C" 

= + C’’ 

and define a discrete-time functional 


(A.l) 


/ N { c^,a = ^\\c 


^IlC 


0\\2 




+ « ^\\Zk - 


2 ^ 


Ilk 


with Mk = AtM when a fixed time-step is considered. This choice of discrete observation 
norm Mk ensures that the discrete-time functional is consistent with respect to the 
continuous-time functional ^. The discrete-time optimal filter is then defined by = 
Alnjargmin _/■„]■ The discrete-time counterpart of (251 for the optimal linear estimator can 


be deduced with a prediction-correction time scheme. We denote by a superscript — the 
prediction computations and by a superscript -I- the correction computations. 

1. Prediction: 



= A„+i|„(X+,0+) 

^n+1 

= S+ 

pxx- 

^n+1 

= (d.A„H.i|„)P- + 


+ (dxA„_|_i|„ 


-f (d^Ayj^xIn, 


+ (deA„+i|„ 

px^- 

-^n+1 

= (dxA„+i|„)Pf + 

-pOO— 

-^n+1 

+ 

II 


(dx 


+ (deA„+i|„)P®®+ 


(A.2a) 


2. Correction: 


n+1 


(P 




r^x 


= (^(dH„+i)TM„+i(dH„+i) 

= P^^+(d,H„+i)TM„+i 
= (Pf++)T(d,H„+i)^M„+i 
= ^n+l + Gl^+l(-^n+l ~ Hn-l-l("^n-|-l)) 
= k+1 + ^n+li^n+1 — 


(A.2b) 


This algorithm can be interpreted as prediction-correction time discretization of (251. 


Another way to propose an efficient estimator is to replace at the time-discrete level 
the tangent operators by a finite difference interpolation based on the computation of the 
original operator on numerous sampling points. This is the case for Ensemble Kalman 
Filter (EnKF) [55] or the Unscented Kalman Filter (UKF) [55j. In this article we focus 
on the second one which is a discrete-time estimator based on sampling particles called 
sigma-points helping to replace the tangent computations. Let us introduce the so-called 
unitary sampling points jl*! and weight Ui with the following rules 


5:a./w=o, 


(A.3) 


2=1 


2=1 
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so that, at each time step, the sigma-points can be generated around the estimated values 
based on the covariance estimation by 


^n+l 

/j[d+ 

V ’^ri+l , 


^;!+i 

Cl 


p+ r[d 


then we compute the prediction 



= A 


n+l|n 
|[^] + 


(aW+,0|:i+), = 


n+1 ’ 

9 - — n-dW” — d+ 


and the corresponding observations 


= H(X 


n-l-ll 


^n+1 — 


2=1 


The gain is then defined by 

= EL 


T)XZ 

-^n+1 


- X-„ 


1 ~ ^^n+l)(^n+l ~ ^n+lJ 

p / r\\i\~ r\— \ / ^\i\~ 


K^+i = ELi - Ci)(C”i - 


pzz 
^n+1 

px _ /'■pxz 

'^n+1 — V-t^n+1 

'^n+1 


= ELi«.(C”i-Cii 

) • (Pf+i)-\ 

= (P^"+i) • (Pf+i)-' 


xzi:!,-- z-^,)T + M, 


r-1 

-‘n+1’ 


(A.4a) 


(A.4b) 


(A.4c) 


(A.4d) 


so that, with the covariance 


p 

'^n+l ~ T! 

i=l 




^n+1 I p0z 
V^n+1 




we still have 


Xn+l — X„_|_i + G^_^_l{Zn+l Z^_^_X, 


Cl = 


n+1 

"P^+i = "P-^, - 'P 


G^ + l(Z„+l - 


n+1 


xz 

n+1 


(P 


ZZ ^ —1 I'tjpxz ^ 'JJ 


n+lJ 


y -L n+i; 


(A.4e) 


This filter has the great advantage to be much easier to compute numerically since, 
contrary to EKF, the tangent operators computations are no longer required. 
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Appendix B. Reduced-Order filters: time discrete RoEKF and RoUKF 

The reduced-order filtering concept presented in Section [T3| can be directly applied 
to the time and space discretized versions of the equations, leading to a time-discrete 
RoEKF: 

1. Prediction: 


2. Correction: 


b V— 
■^n+1 

Ln-t-1 


~ (d ^n+l\rt^^n 


U„+i = \Jn + L^_|_]^(dH„+i)'''M„+i(dH„+i)L„+i 

Gn-i-1 = L„+iU„_|_]^L^_|_]^(dH„_|_i)'''M„+i 

+ Gn+i{Zn+i - 


(B.la) 


(B.lb) 


This algorithm has been applied in [43] for parameter identification by reducing the 
uncertainty space to the parameter space. Indeed, the extension L is initially decomposed 
into 


L(0) = 




Using the fact that the parameters dynamics is null, we can easily proved that 


Vt >0, if = 1. 


(B.2) 


Then, by direct computation, we have the continuous-time formulation of the RoEKF in 
a parameter identification context 


' X = A{X,e)+lie, A(0) = Ao, 

§ = U-iFT(dxH)TM(Z - H(X)), 0(0) = 0o, 

R = (d_,cA)L^ + deA, R(0) = 0, 

U = L^T(d^H)TM(dxH)L^, U(0) = U,. 


In (B.3), we see that the hlter correction appears in the parameters dynamics. 


Then, 

which 


the correction on the parameters is reverberated to the state by the mean of U 
can be interpreted from its proper dynamics as the sensitivity of the model with respect 
to the parameters. Therefore even for a reduced order strategy on the parameters the 
global filter corrects the two components of the model, namely the parameters but also 
the state. 

The same principles can be applied after the time discretization of the model to get 


1. Prediction: 


A-+1 

= A„+i|„(A4,0+) 


^n+l 

= d+ 

(B.4a) 

14+1 

= (dxA„_|_i|„)L„ -I- deA„_|_i|„ 
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2. Correction: 


U„+i = U„ + (L-+i)T(d,H„+OTM„+i(d^H„+i)L-+i 

^n+l — ^n+1 + ^i+l(^n+l ~ ^n+l) 

^n+1 — ^n+1 

n+li^n+l)'^ i^xiin+l)'^^n+liZn+l — H„+i (X„_|_ )) 


(B.4b) 


starting from the same initial conditions. 

We have presented the Reduced Order EKF and have recalled its convergence. Now 
a legitimate question is to know if a reduced order strategy can be applied to other types 
of approximation of the optimal filtering approach. The answer was given for the UKF 
filter in [42] where a reduced order version (RoUKF) was derived. Moreover it was also 
applied to a specihc case where the initial uncertainty is reduced to the parametric space. 
In the general context - i. e. without particularizing the state and parameter dependency 
- the RoUKF can be formulated as follow. 

Given an adequate sampling rule, we store the corresponding weights (at) in the diag¬ 
onal matrix Da and precompute specific unitary simplex sigma-points {P)i<i<p (i.e. with 
zero mean and unit covariance) with p = d-\- \ since the reduced space corresponds here 
to the parametric space. The algorithm consists of the following three steps computed 
recursively: 

1. Sampling: 


rG„ = 

=l+ + L-■C'^/^ 1<*<P 

= 0 ++L^C^/^ 1 < 1 <P 

2. Prediction: 

'I'Vi =A„+i|„(lW+,0t'+), 1<Z<P 

€+1 1<1<P 

^^n+l = Sr=l 

3. Correction: 


■^n+1 

= [X^:l-]Da[I^*^T 


jO 

■^n+1 



^n+1 



^n+l 



Bn-l-l 

= [zW-]Z7„[/W]t 


U„+l 



A++1 


Zn+l) 

^n+l 

= ^n+l + n+l^n+l^n+liZn+l - 

~ Zn+l) 
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(B.5a) 


(B.5b) 


(B.5c) 






where is the matrix concatenating the (/*) vectors side by side, and similarly for 
other vectors. 

The analysis of this estimator can also be found in |32] and relies again on the fact 
that the linearized error satisfies exactly the same dynamics as the linearized error of the 
RoEKF. Therefore the demonstration by linearization is directly obtained. The combined 
practical simplicity and efficiency of such parameter estimator in comparison to other 
adaptive observers have made this approach popular for real case parameter estimation 
problems 
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